Installation Instructions

Download and install miniconda: https://conda.io/miniconda.html

Make sure you are using the conda-forge channel:

$ conda config --add channels conda-forge
$ conda update --yes conda python

Install gsshapy:

$ conda create -n gssha python=2
$ source activate gssha
(gssha)$ conda install --yes gsshapy pynio jupyter

Install GSSHA: http://www.gsshawiki.com/GSSHA_Download

This will NOT work on Windows.

In [17]:
import os
from datetime import datetime, timedelta

from gsshapy.modeling import GSSHAFramework
from gsshapy.grid import HRRRtoGSSHA
from gsshapy.grid.hrrr_to_gssha import download_hrrr_for_gssha
from gsshapy.lib import db_tools as dbt

import pangaea as pa

Setup environment:


In [2]:
# assuming notebook is run from examples folder
# DONT FORGET dos2unix or unix2dos
base_dir = '/Users/rdchlads/GSSHA_INPUT/'
gssha_model_name = '2017_08_16_270m'
gssha_model_directory = os.path.join(base_dir, gssha_model_name)
hrrr_output_directory = os.path.join(gssha_model_directory, 'hrrr_data')
try:
    os.mkdir(hrrr_output_directory)
except OSError:
    pass

Get GSSHA model bounds:


In [3]:
# load in GSSHA model files
project_manager, db_sessionmaker = \
    dbt.get_project_session(gssha_model_name,
                            gssha_model_directory)

db_session = db_sessionmaker()
project_manager.read(directory=gssha_model_directory,
                     filename="{0}.prj".format(gssha_model_name),
                     session=db_session)

gssha_grid = project_manager.getGrid()
# reproject GSSHA grid and get bounds
min_x, max_x, min_y, max_y = gssha_grid.bounds(as_geographic=True)
min_x, max_x, min_y, max_y


Out[3]:
(-123.38788486927815,
 -122.99241054261093,
 38.99598545548028,
 39.399585966942034)

Download HRRR Data:


In [4]:
downloaded_files = download_hrrr_for_gssha(main_directory=hrrr_output_directory,
                                           forecast_start_date_string='20170913',
                                           forecast_start_hour_string='00',
                                           leftlon=min_x, 
                                           rightlon=max_x,
                                           toplat=max_y,
                                           bottomlat=min_y)
downloaded_files


Out[4]:
['/Users/rdchlads/GSSHA_INPUT/2017_08_16_270m/hrrr_data/20170913/hrrr.t00z.wrfsfcf00.grib2',
 '/Users/rdchlads/GSSHA_INPUT/2017_08_16_270m/hrrr_data/20170913/hrrr.t00z.wrfsfcf01.grib2',
 '/Users/rdchlads/GSSHA_INPUT/2017_08_16_270m/hrrr_data/20170913/hrrr.t00z.wrfsfcf02.grib2',
 '/Users/rdchlads/GSSHA_INPUT/2017_08_16_270m/hrrr_data/20170913/hrrr.t00z.wrfsfcf03.grib2',
 '/Users/rdchlads/GSSHA_INPUT/2017_08_16_270m/hrrr_data/20170913/hrrr.t00z.wrfsfcf04.grib2',
 '/Users/rdchlads/GSSHA_INPUT/2017_08_16_270m/hrrr_data/20170913/hrrr.t00z.wrfsfcf05.grib2',
 '/Users/rdchlads/GSSHA_INPUT/2017_08_16_270m/hrrr_data/20170913/hrrr.t00z.wrfsfcf06.grib2',
 '/Users/rdchlads/GSSHA_INPUT/2017_08_16_270m/hrrr_data/20170913/hrrr.t00z.wrfsfcf07.grib2',
 '/Users/rdchlads/GSSHA_INPUT/2017_08_16_270m/hrrr_data/20170913/hrrr.t00z.wrfsfcf08.grib2',
 '/Users/rdchlads/GSSHA_INPUT/2017_08_16_270m/hrrr_data/20170913/hrrr.t00z.wrfsfcf09.grib2',
 '/Users/rdchlads/GSSHA_INPUT/2017_08_16_270m/hrrr_data/20170913/hrrr.t00z.wrfsfcf10.grib2',
 '/Users/rdchlads/GSSHA_INPUT/2017_08_16_270m/hrrr_data/20170913/hrrr.t00z.wrfsfcf11.grib2',
 '/Users/rdchlads/GSSHA_INPUT/2017_08_16_270m/hrrr_data/20170913/hrrr.t00z.wrfsfcf12.grib2',
 '/Users/rdchlads/GSSHA_INPUT/2017_08_16_270m/hrrr_data/20170913/hrrr.t00z.wrfsfcf13.grib2',
 '/Users/rdchlads/GSSHA_INPUT/2017_08_16_270m/hrrr_data/20170913/hrrr.t00z.wrfsfcf14.grib2',
 '/Users/rdchlads/GSSHA_INPUT/2017_08_16_270m/hrrr_data/20170913/hrrr.t00z.wrfsfcf15.grib2',
 '/Users/rdchlads/GSSHA_INPUT/2017_08_16_270m/hrrr_data/20170913/hrrr.t00z.wrfsfcf16.grib2',
 '/Users/rdchlads/GSSHA_INPUT/2017_08_16_270m/hrrr_data/20170913/hrrr.t00z.wrfsfcf17.grib2',
 '/Users/rdchlads/GSSHA_INPUT/2017_08_16_270m/hrrr_data/20170913/hrrr.t00z.wrfsfcf18.grib2']

Inspect the grid data:


In [19]:
with pa.open_mfdataset(downloaded_files,
                       lat_var='gridlat_0',
                       lon_var='gridlon_0',
                       time_var='time',
                       lat_dim='ygrid_0',
                       lon_dim='xgrid_0',
                       time_dim='time',
                       loader='hrrr') as hrd:
    print(hrd)
    print(hrd.PRATE_P0_L1_GLC0)


<xarray.Dataset>
Dimensions:            (time: 19, xgrid_0: 14, ygrid_0: 17)
Coordinates:
    gridlat_0          (ygrid_0, xgrid_0) float32 38.9327 38.9401 38.9476 ...
    gridlon_0          (ygrid_0, xgrid_0) float32 -123.335 -123.301 -123.268 ...
  * time               (time) datetime64[ns] 2017-09-13 2017-09-13T01:00:00 ...
Dimensions without coordinates: xgrid_0, ygrid_0
Data variables:
    TMP_P0_L1_GLC0     (time, ygrid_0, xgrid_0) float64 304.7 304.5 301.9 ...
    RH_P0_L103_GLC0    (time, ygrid_0, xgrid_0) float64 48.6 49.6 54.7 52.8 ...
    PRES_P0_L1_GLC0    (time, ygrid_0, xgrid_0) float64 9.636e+04 9.609e+04 ...
    TMP_P0_L103_GLC0   (time, ygrid_0, xgrid_0) float64 299.7 298.9 297.3 ...
    gridrot_0          (time, ygrid_0, xgrid_0) float32 -0.28069 -0.280328 ...
    TCDC_P0_L10_GLC0   (time, ygrid_0, xgrid_0) float64 5.5 3.625 2.5 2.5 ...
    PRATE_P0_L1_GLC0   (time, ygrid_0, xgrid_0) float64 0.0 0.0 0.0 0.0 0.0 ...
    UGRD_P0_L103_GLC0  (time, ygrid_0, xgrid_0) float64 2.501 2.876 4.001 ...
    VGRD_P0_L103_GLC0  (time, ygrid_0, xgrid_0) float64 0.3044 -0.6331 ...
    DSWRF_P0_L1_GLC0   (time, ygrid_0, xgrid_0) float64 401.8 401.4 405.7 ...
<xarray.DataArray 'PRATE_P0_L1_GLC0' (time: 19, ygrid_0: 17, xgrid_0: 14)>
dask.array<concatenate, shape=(19, 17, 14), dtype=float64, chunksize=(1, 17, 14)>
Coordinates:
    gridlat_0  (ygrid_0, xgrid_0) float32 38.9327 38.9401 38.9476 38.955 ...
    gridlon_0  (ygrid_0, xgrid_0) float32 -123.335 -123.301 -123.268 ...
  * time       (time) datetime64[ns] 2017-09-13 2017-09-13T01:00:00 ...
Dimensions without coordinates: ygrid_0, xgrid_0
Attributes:
    production_status:                              Operational products
    center:                                         US National Weather Servi...
    forecast_time_units:                            hours
    level:                                          [ 0.]
    long_name:                                      Precipitation rate
    parameter_template_discipline_category_number:  [0 0 1 7]
    initial_time:                                   09/13/2017 (00:00)
    grid_type:                                      Lambert Conformal can be ...
    units:                                          kg m-2 s-1
    level_type:                                     Ground or water surface
    parameter_discipline_and_category:              Meteorological products, ...

Map the variable in the GRIB files to the conversion function:


In [7]:
hrrr_forecast_dir = os.path.dirname(downloaded_files[0])
data_var_map_array = [
   ['precipitation_rate', 'PRATE_P0_L1_GLC0'],
   ['pressure', 'PRES_P0_L1_GLC0'],
   ['relative_humidity', 'RH_P0_L103_GLC0'],
   ['wind_speed', ['UGRD_P0_L103_GLC0', 'VGRD_P0_L103_GLC0']],
   ['direct_radiation_cc', ['DSWRF_P0_L1_GLC0', 'TCDC_P0_L10_GLC0']],
   ['diffusive_radiation_cc', ['DSWRF_P0_L1_GLC0', 'TCDC_P0_L10_GLC0']],
   ['temperature', 'TMP_P0_L1_GLC0'],
   ['cloud_cover_pc', 'TCDC_P0_L10_GLC0'],
]

Option 1. Convert Data:


In [18]:
h2g = HRRRtoGSSHA(gssha_project_folder=gssha_model_directory,
                  gssha_project_file_name="{0}.prj".format(gssha_model_name),
                  lsm_input_folder_path=hrrr_forecast_dir,
                  lsm_search_card="hrrr.*.grib2")

# hmet
h2g.lsm_data_to_arc_ascii(data_var_map_array)

# gag
out_gage_file = os.path.join(gssha_model_directory,
                             'gage_hrrr.gag')
h2g.lsm_precip_to_gssha_precip_gage(out_gage_file,
                                    lsm_data_var='PRATE_P0_L1_GLC0',
                                    precip_type='RADAR')

Option 2. Convert Data & Run the model:


In [ ]:
grf = GSSHAFramework("gssha",
                     gssha_model_directory,
                     "{0}.prj".format(gssha_model_name),
                     lsm_folder=hrrr_forecast_dir,
                     lsm_data_var_map_array=data_var_map_array,
                     lsm_precip_data_var='PRATE_P0_L1_GLC0',
                     lsm_precip_type='RADAR',
                     lsm_search_card="hrrr.*.grib2",
                     lsm_lat_var='gridlat_0',
                     lsm_lon_var='gridlon_0',
                     lsm_time_var='time',
                     lsm_lat_dim='ygrid_0',
                     lsm_lon_dim='xgrid_0',
                     lsm_time_dim='time',
                     grid_module='hrrr')

gssha_event_directory = gr.run_forecast()
gssha_event_directory

The gssha_event_directory is where the simulation output is stored.


In [ ]: